%% CLEAR ALL 
clear all;
close all;
clc;
rng(1, 'twister');
tic;

%% ADD PATH 
% Determine where your m-file's folder is.
folder = fileparts(which("DGP_AR1.m")); 
% Add that folder plus all subfolders to the path.
addpath(genpath(folder)); 
% Saving path
file_dir = what('Results'); 
spath=file_dir.path;


% Simulation settings
rho=0.8;        % Autoregressive coefficient 0.3; 0.6 and 0.8
sigma1=1;       % Std of first error term
sigma2=1;       % Std of second error term
alpha=10;       % Degree of misspecification
T = 200;        % Sample size
p=1;            % Lag length
if rho<0.8      % Forecast horizon
hor=10;         % if rho=0.3 or 0.6 then hor=10 
elseif rho>=0.8
hor=20;         % if rho>=0.8 then hor=20
end 
jhor=[0:hor];   % Time line
responseV=2;    % Response variable
normalizeV=1;   % Normalized variable 
weight=0.5;     % Weighting for VAR IRFs
n_boot=500;     % Number of bootstraps
n_MC=1000;       % Number of simulations
blocksize1=4;   % Size of block
blocksize2=20; % Size of block
Ttest=100;      % Size of test set (for cross-validation forecast)

% True IRFs
for j=1:hor+1
IRF_True(:,j)=rho^jhor(:,j);
end

% Run Simulation 
for i_MC=1:n_MC
% Generate error term
epsilon1 = sigma1 * randn(T, 1);
epsilon2 = sigma2 * randn(T, 1);

% Initial values
y = zeros(T, 1);

% Simulate the AR(1) process
for t = 1:T
    if t==1
     y(t) = epsilon1(1)+epsilon2(1);
    elseif t>1
     y(t) = rho*y(t-1)+epsilon1(t)+epsilon2(t)+(alpha/sqrt(T))*epsilon2(t-1); 
    end
end

% Data observed
data=[epsilon1 y];

% Estimate IRFs by VAR
[IRF_VAR_org]=VAR_model(data,hor,p,1,0);
IRF_VAR = IRF_VAR_org(responseV,:)/ IRF_VAR_org(normalizeV,1); 
IRF_VAR_result(:,i_MC)=IRF_VAR';


% Estimate IRFs by LP
[IRF_LP_org]=LP_Jorda_model(data,hor,p,1,0);
IRF_LP = IRF_LP_org(responseV,:)/ IRF_LP_org(normalizeV,1); 
IRF_LP_result(:,i_MC)=IRF_LP';


% Estimate IRFs by Pooling
IRF_Pool=weight*IRF_VAR+(1-weight)*IRF_LP;
IRF_Pool_result(:,i_MC)=IRF_Pool';


% Weighting based on cross-validation forecast
[weightfcVAR, weightfcLP]=Weight_CVforecast(data,p,Ttest,hor);

% Estimate IRFs by Pooling
IRF_Pool_fc=weightfcVAR(responseV,:).*IRF_VAR+weightfcLP(responseV,:).*IRF_LP;
IRF_Pool_fc_result(:,i_MC)=IRF_Pool_fc';


% Sequence of block bootstrap (flexible block size)
[IRF_VAR_boot_org,IRF_LP_boot_org,~,~]=VARLP_seqblockbootstrap(data,p,hor+1,n_boot,responseV,normalizeV,0,0,0);
for j=1:hor+1
covariance_varlp_seqblock(:,:,j)=cov(IRF_VAR_boot_org(:,j),IRF_LP_boot_org(:,j));
correlation_varlp_seqblock(:,:,j)=corrcoef(IRF_VAR_boot_org(:,j),IRF_LP_boot_org(:,j));
correlation_varlp_seqblock_sim(:,:,j,i_MC)=correlation_varlp_seqblock(:,:,j);
end

% Sequence of block bootstrap (fix-block size)
[IRF_VAR_boot_org_bs1,IRF_LP_boot_org_bs1,~,~]=VARLP_seqblockbootstrap(data,p,hor+1,n_boot,responseV,normalizeV,0,1,blocksize1);
[IRF_VAR_boot_org_bs2,IRF_LP_boot_org_bs2,~,~]=VARLP_seqblockbootstrap(data,p,hor+1,n_boot,responseV,normalizeV,0,1,blocksize2);
for j=1:hor+1
covariance_varlp_seqblock_bs1(:,:,j)=cov(IRF_VAR_boot_org_bs1(:,j),IRF_LP_boot_org_bs1(:,j));
correlation_varlp_seqblock_bs1(:,:,j)=corrcoef(IRF_VAR_boot_org_bs1(:,j),IRF_LP_boot_org_bs1(:,j));
correlation_varlp_seqblock_sim_bs1(:,:,j,i_MC)=correlation_varlp_seqblock_bs1(:,:,j);

covariance_varlp_seqblock_bs2(:,:,j)=cov(IRF_VAR_boot_org_bs2(:,j),IRF_LP_boot_org_bs2(:,j));
correlation_varlp_seqblock_bs2(:,:,j)=corrcoef(IRF_VAR_boot_org_bs2(:,j),IRF_LP_boot_org_bs2(:,j));
correlation_varlp_seqblock_sim_bs2(:,:,j,i_MC)=correlation_varlp_seqblock_bs2(:,:,j);
end


% Residual based bootstrap
[IRF_VAR_bootrep,IRF_LP_bootrep]=VAR_bootstrap_replacement(data,hor,p,1,n_boot,responseV,normalizeV);
for j=1:hor+1
covariance_varlp_bootrep(:,:,j)=cov(IRF_VAR_bootrep(:,j),IRF_LP_bootrep(:,j));
correlation_varlp_bootrep(:,:,j)=corrcoef(IRF_VAR_bootrep(:,j),IRF_LP_bootrep(:,j));
correlation_varlp_bootrep_sim(:,:,j,i_MC)=correlation_varlp_bootrep(:,:,j);
end

scale_VAR =1/ IRF_VAR_org(normalizeV,1); 
scale_LP =1/ IRF_LP_org(normalizeV,1);
zvalue=1.96;

for j=1:hor+1
% For sequence of block bootstrap (flexible block size)
IRF_VAR_variance_seqblock(:,j)=covariance_varlp_seqblock(1,1,j);
IRF_LP_variance_seqblock(:,j)=covariance_varlp_seqblock(2,2,j);
IRF_VARLP_covariance_seqblock(:,j)=covariance_varlp_seqblock(2,1,j);

% For sequence of block bootstrap (fix-block size)
IRF_VAR_variance_seqblock_bs1(:,j)=covariance_varlp_seqblock_bs1(1,1,j);
IRF_LP_variance_seqblock_bs1(:,j)=covariance_varlp_seqblock_bs1(2,2,j);
IRF_VARLP_covariance_seqblock_bs1(:,j)=covariance_varlp_seqblock_bs1(2,1,j);

IRF_VAR_variance_seqblock_bs2(:,j)=covariance_varlp_seqblock_bs2(1,1,j);
IRF_LP_variance_seqblock_bs2(:,j)=covariance_varlp_seqblock_bs2(2,2,j);
IRF_VARLP_covariance_seqblock_bs2(:,j)=covariance_varlp_seqblock_bs2(2,1,j);

% For residual based bootstrap
IRF_VAR_variance_bootrep(:,j)=covariance_varlp_bootrep(1,1,j);
IRF_LP_variance_bootrep(:,j)=covariance_varlp_bootrep(2,2,j);
IRF_VARLP_covariance_bootrep(:,j)=covariance_varlp_bootrep(2,1,j);


% Variance of Pooled IRFs
IRF_Pool_sd_seqblock(:,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_seqblock(:,j)+(1-weight)^2*scale_LP^2*IRF_LP_variance_seqblock(:,j)...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_seqblock(:,j));

IRF_Pool_sd_seqblock_bs1(:,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_seqblock_bs1(:,j)+(1-weight)^2*scale_LP^2*IRF_LP_variance_seqblock_bs1(:,j)...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_seqblock_bs1(:,j));

IRF_Pool_sd_seqblock_bs2(:,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_seqblock_bs2(:,j)+(1-weight)^2*scale_LP^2*IRF_LP_variance_seqblock_bs2(:,j)...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_seqblock_bs2(:,j));

IRF_Pool_sd_bootrep(:,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_bootrep(:,j)+(1-weight)^2*scale_LP^2*IRF_LP_variance_bootrep(:,j)...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_bootrep(:,j));
end

% Confidence intervals and lengths
% Confidence intervals Sequence Block Bootstrap (flexible block size)
IRF_VAR_UB_seqblock(:,i_MC)=(IRF_VAR+zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_seqblock))';
IRF_VAR_LB_seqblock(:,i_MC)=(IRF_VAR-zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_seqblock))';

IRF_LP_UB_seqblock(:,i_MC)=(IRF_LP+zvalue*sqrt(scale_LP^2*IRF_LP_variance_seqblock))';
IRF_LP_LB_seqblock(:,i_MC)=(IRF_LP-zvalue*sqrt(scale_LP^2*IRF_LP_variance_seqblock))';

IRF_Pool_UB_seqblock(:,i_MC)=(IRF_Pool+zvalue*IRF_Pool_sd_seqblock)';
IRF_Pool_LB_seqblock(:,i_MC)=(IRF_Pool-zvalue*IRF_Pool_sd_seqblock)';

% Confidence intervals Sequence Block Bootstrap (fix block size)
IRF_Pool_UB_seqblock_bs1(:,i_MC)=(IRF_Pool+zvalue*IRF_Pool_sd_seqblock_bs1)';
IRF_Pool_LB_seqblock_bs1(:,i_MC)=(IRF_Pool-zvalue*IRF_Pool_sd_seqblock_bs1)';

IRF_Pool_UB_seqblock_bs2(:,i_MC)=(IRF_Pool+zvalue*IRF_Pool_sd_seqblock_bs2)';
IRF_Pool_LB_seqblock_bs2(:,i_MC)=(IRF_Pool-zvalue*IRF_Pool_sd_seqblock_bs2)';


% Confidence intervals Bootstrap with replacement
IRF_VAR_UB_bootrep(:,i_MC)=(IRF_VAR+zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_bootrep))';
IRF_VAR_LB_bootrep(:,i_MC)=(IRF_VAR-zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_bootrep))';

IRF_LP_UB_bootrep(:,i_MC)=(IRF_LP+zvalue*sqrt(scale_LP^2*IRF_LP_variance_bootrep))';
IRF_LP_LB_bootrep(:,i_MC)=(IRF_LP-zvalue*sqrt(scale_LP^2*IRF_LP_variance_bootrep))';

IRF_Pool_UB_bootrep(:,i_MC)=(IRF_Pool+zvalue*IRF_Pool_sd_bootrep)';
IRF_Pool_LB_bootrep(:,i_MC)=(IRF_Pool-zvalue*IRF_Pool_sd_bootrep)';

% Confidence interval length
IRF_VAR_seqblock_length(:,i_MC)=(2*zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_seqblock))';
IRF_LP_seqblock_length(:,i_MC)=(2*zvalue*sqrt(scale_LP^2*IRF_LP_variance_seqblock))';
IRF_Pool_seqblock_length(:,i_MC)=(2*zvalue*IRF_Pool_sd_seqblock)';
IRF_Pool_seqblock_length_bs1(:,i_MC)=(2*zvalue*IRF_Pool_sd_seqblock_bs1)';
IRF_Pool_seqblock_length_bs2(:,i_MC)=(2*zvalue*IRF_Pool_sd_seqblock_bs2)';


IRF_VAR_bootrep_length(:,i_MC)=(2*zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_bootrep))';
IRF_LP_bootrep_length(:,i_MC)=(2*zvalue*sqrt(scale_LP^2*IRF_LP_variance_bootrep))';
IRF_Pool_bootrep_length(:,i_MC)=(2*zvalue*IRF_Pool_sd_bootrep)';

end % End simulation


% Estimate all results
Run_ARDGP_all_results

% Save results
Run_ARDGP_save_results
